\(~\)
The pigeon_video.mp4 file and priming emails sent to
students are available in the Teaching_resources folder
found on the github page associated with this manuscript. They can be
accessed here: https://github.com/tomkeaney/Biased_pigeons.
‘Hungry forager’ primer: “The video was recorded just after winter, which is typically the most difficult season of the year for pigeons in terms of feeding and survival. About 30% of pigeons do not survive winter; many that do survive experience a significant loss of body mass and must make the most of every feeding opportunity that arises in order to regain condition.”
‘Satiated forager’ primer: “The video was recorded late in the afternoon near Jewell railway station, a favourite spot for commuters to feed the pigeons. Moreland Council estimates that locals dump 5-10kg of bread and seed here, creating a superabundant daily food source. By the end of the day, many of the birds are completely satiated and have little or no motivation to feed.”
\(~\)
library(tidyverse) # tidy style coding
library(brms) # Bayesian models
library(loo) # for information criteria
library(tidybayes) # Bayesian aesthetics
library(MetBrewer) # colours
library(kableExtra) # tables
library(patchwork) # putting plots together
library(DT) # for search- and saveable tables
library(pander) # for simpler tables
library(png) # to load images
library(grid) # to plot images
library(ggdag) # to draw dags
We collated 89 responses across three years, from 82 students. Four students from the 2021 cohort provided multiple responses; we chose to only include their initial response as we cannot be sure that later responses were made while the students were still unaware of the experiment. One additional response did not provide the indicated expectation of the student, so it was removed from the dataset. The cleaned dataset therefore contains 81 responses from 81 students.
De-identified data can be downloaded directly from the table below (click the CSV button), or from the data folder in the github repository associated with this report.
data <-
read_csv("data/pigeon_bias_data.csv") %>%
mutate(Student_ID = as.factor(Student_ID),
Year = as.factor(Year),
Foraging_prop = (Foraging_percentage / 100)) %>%
filter(First_response == "YES",
Foraging_percentage != "NA",
Primer_understood != "NA") %>%
select(Student_ID, Year, First_response, everything())
data_peck <-
data %>%
filter(Peck_rate_2 != "NA") %>%
pivot_longer(cols = Peck_rate_1:Peck_rate_2, names_to = "Trial",
values_to = "Peck_rate")
# Create a function to build HTML searchable tables
my_data_table <- function(df){
datatable(
df, rownames=FALSE,
autoHideNavigation = TRUE,
extensions = c("Scroller", "Buttons"),
options = list(
dom = 'Bfrtip',
deferRender=TRUE,
scrollX=TRUE, scrollY=400,
scrollCollapse=TRUE,
buttons =
list('pageLength', 'colvis', 'csv', list(
extend = 'pdf',
pageSize = 'A4',
orientation = 'landscape',
filename = 'full_dataset')),
pageLength = 81
)
)
}
my_data_table(data)
Column explanations
Student_ID: unique, anonymised identifier for each student
Year: year that the experiment was conducted.
Gender: upon enrollment at The University of Melbourne, students are asked to indicate their title. We identified women as observers that answered “Miss” or “Ms” and Men as those who that answered “Mr”. Those with entirely missing entries were coded as “NA”.
First_response: indicates whether this was the initial response provided by the student.
Bias_treatment: the primer the observer received, where ‘satiated’ indicates that the observers were provided information prior to a trial that suggested pigeons were fell fed, whereas ‘hungry’ indicated that the pigeons were in poor condition and hungry.
Expectation: we asked the observers to indicate whether they thought the pigeons would be hungry or satiated. We included this question to test whether observers were appropriately primed by their bias treatment.
Primer_understood: did the observer’s expectation match the primer they received?
Foraging_percentage: the percentage of pigeons that observers estimated to be foraging over a 15 second period, while observing a large flock.
Peck_rate_1: the number of times a single chosen pigeon pecked the ground over a 15 second period.
Peck_rate_2: the number of times a second chosen pigeon pecked the ground over a 15 second period.
Foraging_prop: proportion of pigeons estimated to be foraging.
\(~\)
\(~\)
\(~\)
We find that 24 of the 81 observers indicated a feeding motivation expectation opposite to that implied by the primer they were allocated.
This suggests that Expectation may be a better predictor
of foraging estimation than allocated Bias_treatment.
We explicitly assess the effect of treatment and expectation by fitting two models:
a model with allocated primer (Bias_treatment) as
the predictor variable
a model with indicated hunger expectation
(Expectation) as the predictor variable
\(~\)
# First let's model the effect of bias treatment on foraging estimation
foraging_model_treatment <- brm(Foraging_prop ~ 0 + Bias_treatment + Year,
data = data, family = Beta,
prior = c(prior(normal(0, 1.5), class = b),
prior(exponential(1), class = phi)),
iter = 6000, warmup = 2000, chains = 4, cores = 4,
control = list(adapt_delta = 0.8, max_treedepth = 10),
seed = 1, file = "fits/foraging_model_treatment")
foraging_model_treatment
## Family: beta
## Links: mu = logit; phi = identity
## Formula: Foraging_prop ~ 0 + Bias_treatment + Year
## Data: data (Number of observations: 81)
## Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 16000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Bias_treatmentHungry -0.43 0.22 -0.86 -0.01 1.00 6351
## Bias_treatmentSatiated -0.40 0.23 -0.86 0.05 1.00 6299
## Year2022 -0.02 0.25 -0.51 0.48 1.00 6471
## Year2023 -0.21 0.25 -0.70 0.29 1.00 6240
## Tail_ESS
## Bias_treatmentHungry 8461
## Bias_treatmentSatiated 8510
## Year2022 8733
## Year2023 9126
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 4.23 0.60 3.13 5.51 1.00 11656 10500
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Table S1. Posterior estimates of the percentage of pigeons foraging, split by the primer studnets were allocated.
new_data <- expand_grid(Bias_treatment = c("Hungry", "Satiated"),
Year = 2023)
new_data %>%
left_join(data %>% group_by(Bias_treatment) %>% summarise(`n students` = n())) %>%
cbind(fitted(foraging_model_treatment, newdata = new_data, summary = T) %>%
as_tibble() %>%
mutate(across(1:4, ~ .x *100),
across(1:4, round, 2))) %>%
rename("Estimated % foraging" = Estimate,
"Bias treatment" = Bias_treatment) %>%
select(-Year) %>%
pander(split.cell = 20, split.table = Inf)
| Bias treatment | n students | Estimated % foraging | Est.Error | Q2.5 | Q97.5 |
|---|---|---|---|---|---|
| Hungry | 44 | 34.84 | 4.32 | 26.7 | 43.68 |
| Satiated | 37 | 35.36 | 4.27 | 27.24 | 44.09 |
# fit the same model, except using participant expectation rather than allocated bias treatment
foraging_model_expectation <- brm(Foraging_prop ~ 0 + Expectation + Year,
data = data, family = Beta,
prior = c(prior(normal(0, 1.5), class = b),
prior(exponential(1), class = phi)),
iter = 6000, warmup = 2000, chains = 4, cores = 4,
control = list(adapt_delta = 0.8, max_treedepth = 10),
seed = 1, file = "fits/foraging_model_expectation")
foraging_model_expectation
## Family: beta
## Links: mu = logit; phi = identity
## Formula: Foraging_prop ~ 0 + Expectation + Year
## Data: data (Number of observations: 81)
## Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 16000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ExpectationHungry -0.29 0.23 -0.75 0.15 1.00 5768 8581
## ExpectationSatiated -0.51 0.21 -0.93 -0.09 1.00 6442 8860
## Year2022 -0.06 0.25 -0.55 0.44 1.00 5845 8238
## Year2023 -0.24 0.26 -0.74 0.27 1.00 5882 7991
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 4.29 0.61 3.19 5.58 1.00 12560 11169
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Table S2. Posterior estimates of the percentage of pigeons foraging, split by the actual expectation of observers.
new_data_2 <- expand_grid(Expectation = c("Hungry", "Satiated"),
Year = 2023)
new_data_2 %>%
left_join(data %>% group_by(Expectation) %>% summarise(`n observers` = n())) %>%
cbind(fitted(foraging_model_expectation, newdata = new_data_2, summary = T) %>%
as_tibble() %>%
mutate(across(1:4, ~ .x *100),
across(1:4, round, 2))) %>%
rename("Estimated % foraging" = Estimate,
"Indicated expectation" = Expectation) %>%
select(-Year) %>%
pander(split.cell = 20, split.table = Inf)
| Indicated expectation | n observers | Estimated % foraging | Est.Error | Q2.5 | Q97.5 |
|---|---|---|---|---|---|
| Hungry | 45 | 36.98 | 4.17 | 29.13 | 45.32 |
| Satiated | 36 | 32.24 | 4.49 | 23.75 | 41.32 |
\(~\)
\(~\)
Get posterior means and difference contrasts
# treatment model
draws_treatment <-
as_draws_df(foraging_model_treatment) %>%
mutate(Hungry = inv_logit_scaled(b_Bias_treatmentHungry + b_Year2023) *100,
Satiated = inv_logit_scaled(b_Bias_treatmentSatiated + b_Year2023)*100,
diff_contrast = (Hungry - Satiated)) %>%
select(Hungry, Satiated, diff_contrast) %>%
pivot_longer(names_to = "Treatment", values_to = "Posterior_estimate", cols = 1:3) %>%
mutate(Predictor = "Allocated primer")
p2 <-
draws_treatment %>%
filter(Treatment != "diff_contrast") %>%
ggplot(aes(x = Treatment, y = Posterior_estimate)) +
stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5) +
scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
coord_flip() +#ylim = c(25, 55)) +
#geom_vline(xintercept = 0, linetype = 2) +
#scale_y_continuous(breaks = c(, 0, 1)) +
xlab("Allocated primer") +
ylab("Estimated % pigeons foraging") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size = 14))
p3 <-
draws_treatment %>%
filter(Treatment == "diff_contrast") %>%
ggplot(aes(y = Posterior_estimate)) +
stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5, scale =0.5) +
scale_fill_manual(values = met.brewer("Hiroshige")[4]) +
coord_flip(ylim = c(-20, 20)) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
#scale_y_continuous(breaks = c(, 0, 1)) +
xlab(NULL) +
ylab("Hungry - Satiated difference\ncontrast (% points)") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
text = element_text(size = 14))
# expectation model
draws_expectation <-
as_draws_df(foraging_model_expectation) %>%
mutate(Hungry = inv_logit_scaled(b_ExpectationHungry + b_Year2023) *100,
Satiated = inv_logit_scaled(b_ExpectationSatiated + b_Year2023)*100,
diff_contrast = (Hungry - Satiated)) %>%
select(Hungry, Satiated, diff_contrast) %>%
pivot_longer(names_to = "Treatment", values_to = "Posterior_estimate", cols = 1:3) %>%
mutate(Predictor = "Indicated expectation")
p4 <-
draws_expectation %>%
filter(Treatment != "diff_contrast") %>%
ggplot(aes(x = Treatment, y = Posterior_estimate)) +
stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5) +
scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
coord_flip() + #ylim = c(25, 55)) +
#geom_vline(xintercept = 0, linetype = 2) +
#scale_y_continuous(breaks = c(, 0, 1)) +
xlab("Indicated expectation") +
ylab("Estimated % pigeons foraging") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size = 14))
p5 <-
draws_expectation %>%
filter(Treatment == "diff_contrast") %>%
ggplot(aes(y = Posterior_estimate)) +
stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5, scale =0.5) +
scale_fill_manual(values = met.brewer("Hiroshige")[4]) +
coord_flip(ylim = c(-20, 20)) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
#scale_y_continuous(breaks = c(, 0, 1)) +
xlab(NULL) +
ylab("Hungry - Satiated difference\ncontrast (% points)") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
text = element_text(size = 14))
\(~\)
\(~\)
We asked observers to count the number of pecks of the ground that a selected pigeon made over a 1 minute period. We use the number of pecks that occur as a measure of feeding rate.
We first estimated a baseline peck rate by observing 35 pigeons. To select the pigeons we observed, we split a still image of the foraging video (taken at time zero) into a 43 x 21 cell grid. From the 122 cells that contained pigeons, 40 were chosen by random number generation (see code chunk below). In the event that multiple pigeons were present in the cell, we selected the most prominent to observe. Five observations were discarded - three due to overlap of the same pigeon between cells that were selected by the random number generator, and two more as the pigeons left the field of view during the video and could no longer be tracked.
# curly brackets run all lines included within them
{set.seed(1) # so that sample produces a reproducible sequence
sample(1:122, 40, replace = FALSE)}
## [1] 121 68 39 1 34 87 43 14 82 59 51 97 85 21 106 54 74 7 73
## [20] 79 110 37 89 101 118 100 44 103 33 84 35 70 108 42 38 20 28 117
## [39] 96 91
The selected pigeons for baseline observation are shown in the image below
img <- readPNG("pigeon_selection.png")
grid.raster(img)
\(~\)
\(~\)
Load in the data
baseline_data <- read_csv("data/baseline_peck_data.csv") %>%
select(1:5) %>% # remove the comments column
pivot_longer(cols = 4:5, names_to = "Observation", values_to = "Peck_rate") %>%
mutate(ID = as.factor(ID),
Observation = str_remove(Observation, "Peck_count_")) %>%
rename(Pigeon_ID = ID) %>%
filter(!is.na(Peck_rate))
my_data_table <- function(df){
datatable(
df, rownames=FALSE,
autoHideNavigation = TRUE,
extensions = c("Scroller", "Buttons"),
options = list(
dom = 'Bfrtip',
deferRender=TRUE,
scrollX=TRUE, scrollY=400,
scrollCollapse=TRUE,
buttons =
list('pageLength', 'colvis', 'csv', list(
extend = 'pdf',
pageSize = 'A4',
orientation = 'landscape',
filename = 'baseline_dataset')),
pageLength = 78
)
)
}
my_data_table(baseline_data)
X and Y represent grid
coordinates.
Pigeon_ID identifies a specific pigeon
Observation indicates whether this was the first or
second scoring for a single pigeon. We scored each pigeon twice as
distant pigeons were difficult to observe and to ensure that the correct
pigeon was tracked throughout the minute of observation.
Peck_rate is the number of times the ground was
pecked over a minute of observation.
\(~\)
Fit a simple model to estimate median peck rate
baseline_peck_model_zi <-
brm(Peck_rate ~ 1 + (1|Pigeon_ID),
family = zero_inflated_negbinomial(), data = baseline_data,
prior = c( prior(normal(0, 1.5), class = Intercept),
prior(exponential(1), class = sd),
prior(exponential(1), class = shape),
prior(exponential(1), class = zi)),
chains = 4, cores = 4, warmup = 2000, iter = 6000,
file = "fits/baseline_peck_model")
# wrangle the output
baseline_peck_predictions <-
baseline_peck_model_zi %>%
as_draws_df() %>%
mutate(Baseline_estimate = exp(b_Intercept),
peck_rate_sd = exp(sd_Pigeon_ID__Intercept)) %>%
select(Baseline_estimate, peck_rate_sd)
baseline_data %>%
distinct(Pigeon_ID) %>%
summarise(`n pigeons observed` = length(Pigeon_ID)) %>%
bind_cols(
fitted(baseline_peck_model_zi, summary = T, re_formula = NA) %>%
as_tibble() %>%
distinct(Estimate, .keep_all = T) %>%
mutate(across(1:4, round, 2)) %>%
rename(`Baseline median peck rate / per min` = Estimate)) %>%
pander()
| n pigeons observed | Baseline median peck rate / per min | Est.Error | Q2.5 | Q97.5 |
|---|---|---|---|---|
| 35 | 1.32 | 0.51 | 0.54 | 2.53 |
\(~\)
\(~\)
\(~\)
Once again, we hypothesised that attention paid to and/or comprehension of the primer statement had a large effect on observers’ perception of pigeon foraging.
Lets again fit the two models:
a model with allocated primer (Bias_treatment) as
the predictor variable
a model with indicated hunger expectation
(Expectation) as the predictor variable
\(~\)
peck_model_treatment <-
brm(Peck_rate ~ 0 + Bias_treatment + Year + (1|Student_ID),
data = data_peck, family = negbinomial,
prior = c(prior(normal(0, 1.5), class = b),
prior(exponential(1), class = sd)),
iter = 6000, warmup = 2000, chains = 4, cores = 4,
control = list(adapt_delta = 0.95, max_treedepth = 12),
seed = 1, file = "fits/peck_model_treatment")
peck_model_treatment
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: Peck_rate ~ 0 + Bias_treatment + Year + (1 | Student_ID)
## Data: data_peck (Number of observations: 162)
## Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 16000
##
## Group-Level Effects:
## ~Student_ID (Number of levels: 81)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.22 0.15 0.01 0.55 1.00 4301 6867
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Bias_treatmentHungry 2.30 0.20 1.92 2.69 1.00 9114
## Bias_treatmentSatiated 2.13 0.22 1.70 2.57 1.00 9172
## Year2022 0.03 0.24 -0.45 0.50 1.00 9925
## Year2023 0.01 0.24 -0.47 0.48 1.00 9503
## Tail_ESS
## Bias_treatmentHungry 11051
## Bias_treatmentSatiated 11147
## Year2022 11314
## Year2023 11012
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 0.88 0.12 0.68 1.13 1.00 12356 10493
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Table S3. The estimated peck rate of foraging pigeons, split by the primer observers were allocated.
new_data %>%
bind_cols(fitted(peck_model_treatment, newdata = new_data, summary = T, re_formula = NA) %>%
as_tibble() %>%
mutate(across(1:4, round, 2))) %>%
left_join(data_peck %>% group_by(Bias_treatment) %>%
distinct(Student_ID) %>% summarise(`n pigeons observed` = n())) %>%
rename("Estimated peck rate" = Estimate,
"Bias treatment" = Bias_treatment) %>%
bind_rows(fitted(baseline_peck_model_zi, summary = T, re_formula = NA) %>%
as_tibble() %>%
distinct(Estimate, .keep_all = T) %>%
mutate(across(1:4, round, 2)) %>%
rename("Estimated peck rate" = Estimate) %>%
mutate(`Bias treatment` = "Baseline") %>%
bind_cols(baseline_data %>%
distinct(Pigeon_ID) %>%
summarise(`n pigeons observed` = length(Pigeon_ID)))) %>%
select(`Bias treatment`, `n pigeons observed`, `Estimated peck rate`, Est.Error, Q2.5, Q97.5) %>%
pander(split.cell = 20, split.table = Inf)
| Bias treatment | n pigeons observed | Estimated peck rate | Est.Error | Q2.5 | Q97.5 |
|---|---|---|---|---|---|
| Hungry | 44 | 10.25 | 1.91 | 7.02 | 14.48 |
| Satiated | 37 | 8.64 | 1.62 | 5.94 | 12.26 |
| Baseline | 35 | 1.32 | 0.51 | 0.54 | 2.53 |
# fit the same model, except using participant expectation rather than allocated bias treatment
peck_model_expectation <- brm(Peck_rate ~ 0 + Expectation + Year + (1|Student_ID),
data = data_peck, family = negbinomial,
prior = c(prior(normal(0, 1.5), class = b),
prior(exponential(1), class = sd)),
iter = 6000, warmup = 2000, chains = 4, cores = 4,
control = list(adapt_delta = 0.95, max_treedepth = 12),
seed = 1, file = "fits/peck_model_expectation")
peck_model_expectation
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: Peck_rate ~ 0 + Expectation + Year + (1 | Student_ID)
## Data: data_peck (Number of observations: 162)
## Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
## total post-warmup draws = 16000
##
## Group-Level Effects:
## ~Student_ID (Number of levels: 81)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.19 0.14 0.01 0.50 1.00 4661 7082
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ExpectationHungry 2.45 0.20 2.07 2.85 1.00 10109 10373
## ExpectationSatiated 1.99 0.21 1.60 2.40 1.00 10768 11307
## Year2022 -0.04 0.23 -0.50 0.42 1.00 10833 11292
## Year2023 -0.03 0.24 -0.50 0.43 1.00 10141 10468
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 0.90 0.12 0.70 1.16 1.00 14648 10643
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Table S4. The estimated peck rate of foraging pigeons, split by the indicated expectation of the observers.
new_data_2 %>%
bind_cols(fitted(peck_model_expectation, newdata = new_data_2, summary = T, re_formula = NA) %>%
as_tibble() %>%
mutate(across(1:4, round, 2))) %>%
left_join(data_peck %>% group_by(Expectation) %>%
distinct(Student_ID) %>% summarise(`n pigeons observed` = n())) %>%
rename("Estimated peck rate" = Estimate,
"Indicated expectation" = Expectation) %>%
bind_rows(fitted(baseline_peck_model_zi, summary = T, re_formula = NA) %>%
as_tibble() %>%
distinct(Estimate, .keep_all = T) %>%
mutate(across(1:4, round, 2)) %>%
rename("Estimated peck rate" = Estimate) %>%
mutate(`Indicated expectation` = "Baseline") %>%
bind_cols(baseline_data %>%
distinct(Pigeon_ID) %>%
summarise(`n pigeons observed` = length(Pigeon_ID)))) %>%
select(`Indicated expectation`, `n pigeons observed`, `Estimated peck rate`, Est.Error, Q2.5, Q97.5) %>%
pander(split.cell = 20, split.table = Inf)
| Indicated expectation | n pigeons observed | Estimated peck rate | Est.Error | Q2.5 | Q97.5 |
|---|---|---|---|---|---|
| Hungry | 45 | 11.44 | 2.01 | 8.05 | 15.94 |
| Satiated | 36 | 7.23 | 1.36 | 4.94 | 10.25 |
| Baseline | 35 | 1.32 | 0.51 | 0.54 | 2.53 |
\(~\)
\(~\)
Get posterior means and difference contrasts
# treatment model
peck_draws_treatment <-
as_draws_df(peck_model_treatment) %>%
mutate(Hungry = exp(b_Bias_treatmentHungry + b_Year2023),
Satiated = exp(b_Bias_treatmentSatiated + b_Year2023),
diff_contrast = (Hungry - Satiated)) %>%
select(Hungry, Satiated, diff_contrast) %>%
bind_cols(baseline_peck_predictions %>% select(Baseline_estimate)) %>%
pivot_longer(names_to = "Treatment", values_to = "Posterior_estimate", cols = 1:3) %>%
mutate(Predictor = "Allocated primer")
p6 <-
peck_draws_treatment %>%
filter(Treatment != "diff_contrast") %>%
ggplot(aes(x = Treatment, y = Posterior_estimate)) +
stat_slab(aes(y = Baseline_estimate),
linetype = 2, linewidth = 0.8, slab_fill = "white",
colour = "black") +
stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5) +
scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
coord_flip(ylim = c(0, 20)) +
#geom_vline(xintercept = 0, linetype = 2) +
#scale_y_continuous(breaks = c(, 0, 1)) +
xlab("Allocated primer") +
ylab("Estimated pecks per min") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size = 14))
p7 <-
peck_draws_treatment %>%
filter(Treatment == "diff_contrast") %>%
ggplot(aes(y = Posterior_estimate)) +
stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5, scale =0.5) +
scale_fill_manual(values = met.brewer("Hiroshige")[4]) +
coord_flip(ylim = c(-5, 12)) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
#scale_y_continuous(breaks = c(, 0, 1)) +
xlab(NULL) +
ylab("Hungry - Satiated difference\ncontrast (pecks per min)") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
text = element_text(size = 14))
# expectation model
peck_draws_expectation <-
as_draws_df(peck_model_expectation) %>%
mutate(Hungry = exp(b_ExpectationHungry + b_Year2023),
Satiated = exp(b_ExpectationSatiated + b_Year2023),
diff_contrast = (Hungry - Satiated)) %>%
select(Hungry, Satiated, diff_contrast) %>%
bind_cols(baseline_peck_predictions %>% select(Baseline_estimate)) %>%
pivot_longer(names_to = "Treatment", values_to = "Posterior_estimate", cols = 1:3) %>%
mutate(Predictor = "Allocated primer")
p8 <-
peck_draws_expectation %>%
filter(Treatment != "diff_contrast") %>%
ggplot(aes(x = Treatment, y = Posterior_estimate)) +
stat_slab(aes(y = Baseline_estimate),
linetype = 2, linewidth = 0.8, slab_fill = "white",
colour = "black") +
stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5) +
scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
coord_flip(ylim = c(0, 20)) +
#geom_vline(xintercept = 0, linetype = 2) +
#scale_y_continuous(breaks = c(, 0, 1)) +
xlab("Indicated expectation") +
ylab("Estimated pecks per min") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size = 14))
p9 <-
peck_draws_expectation %>%
filter(Treatment == "diff_contrast") %>%
ggplot(aes(y = Posterior_estimate)) +
stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.5, scale =0.5) +
scale_fill_manual(values = met.brewer("Hiroshige")[4]) +
coord_flip(ylim = c(-5, 12)) +
geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
#scale_y_continuous(breaks = c(, 0, 1)) +
xlab(NULL) +
ylab("Hungry - Satiated difference\ncontrast (pecks per min)") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
text = element_text(size = 14))
\(~\)
Table S5. The degree to which each group of observer’s overestimates feeding rate (number of ground pecks per minute)
baseline_peck_predictions %>% select(Baseline_estimate) %>% bind_cols(
as_draws_df(peck_model_treatment) %>%
mutate(Hungry = exp(b_Bias_treatmentHungry + b_Year2023),
Satiated = exp(b_Bias_treatmentSatiated + b_Year2023)) %>%
select(Hungry, Satiated)) %>%
mutate(`Bias treatment Satiated / Baseline` = Satiated / Baseline_estimate,
`Bias treatment Hungry / Baseline` = Hungry / Baseline_estimate) %>%
select(contains("Bias")) %>%
pivot_longer(cols = everything(), values_to = "estimate", names_to = "Stat") %>%
group_by(Stat) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE), .cores = 4) %>%
select(-variable) %>%
bind_rows(
baseline_peck_predictions %>% select(Baseline_estimate) %>% bind_cols(
as_draws_df(peck_model_expectation) %>%
mutate(Hungry = exp(b_ExpectationHungry + b_Year2023),
Satiated = exp(b_ExpectationSatiated + b_Year2023)) %>%
select(Hungry, Satiated)) %>%
mutate(`Expectation Satiated / Baseline` = Satiated / Baseline_estimate,
`Expectation Hungry / Baseline` = Hungry / Baseline_estimate) %>%
select(contains("Expectation")) %>%
pivot_longer(cols = everything(), values_to = "estimate", names_to = "Stat") %>%
group_by(Stat) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE), .cores = 4) %>%
select(-variable)
) %>%
pander()
| Stat | median | sd | 2.5% | 97.5% |
|---|---|---|---|---|
| Bias treatment Hungry / Baseline | 7.835 | 4.344 | 3.532 | 19.61 |
| Bias treatment Satiated / Baseline | 6.582 | 3.649 | 2.967 | 16.49 |
| Expectation Hungry / Baseline | 8.74 | 4.825 | 4.013 | 22.16 |
| Expectation Satiated / Baseline | 5.496 | 3.099 | 2.466 | 13.92 |
Option 1
(p2 + p3) / (p4 + p5) / (p6 + p7) / (p8 + p9) +
plot_annotation(tag_levels = 'a')
Figure 1. Posterior mean estimates and difference contrasts for observer estimated group foraging percentage and individual feeding rates. The coloured area is the posterior distribution and the white point is the mean estimate with associated 67% and 95% credible intervals. The distributions shown with dashed lines are the posterior for baseline peck rate, estimated from a random sample of pigeons appearing in the video.
sessionInfo() %>% pander
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
locale: LC_COLLATE=English_Australia.utf8, LC_CTYPE=English_Australia.utf8, LC_MONETARY=English_Australia.utf8, LC_NUMERIC=C and LC_TIME=English_Australia.utf8
attached base packages: grid, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: ggdag(v.0.2.10), png(v.0.1-8), pander(v.0.6.5), DT(v.0.28), patchwork(v.1.1.2), kableExtra(v.1.3.4), MetBrewer(v.0.2.0), tidybayes(v.3.0.4), loo(v.2.6.0), brms(v.2.19.0), Rcpp(v.1.0.11), lubridate(v.1.9.2), forcats(v.1.0.0), stringr(v.1.5.0), dplyr(v.1.1.2), purrr(v.1.0.1), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.2) and tidyverse(v.2.0.0)
loaded via a namespace (and not attached): tensorA(v.0.36.2), rstudioapi(v.0.15.0), jsonlite(v.1.8.7), magrittr(v.2.0.3), farver(v.2.1.1), rmarkdown(v.2.23), vctrs(v.0.6.3), base64enc(v.0.1-3), webshot(v.0.5.5), htmltools(v.0.5.5), distributional(v.0.3.2), curl(v.5.0.1), sass(v.0.4.7), StanHeaders(v.2.26.27), bslib(v.0.5.0), htmlwidgets(v.1.6.2), plyr(v.1.8.8), zoo(v.1.8-12), cachem(v.1.0.8), igraph(v.1.5.0.1), mime(v.0.12), lifecycle(v.1.0.3), pkgconfig(v.2.0.3), colourpicker(v.1.2.0), Matrix(v.1.5-4.1), R6(v.2.5.1), fastmap(v.1.1.1), shiny(v.1.7.4.1), digest(v.0.6.33), colorspace(v.2.1-0), ps(v.1.7.5), crosstalk(v.1.2.0), labeling(v.0.4.2), fansi(v.1.0.4), timechange(v.0.2.0), httr(v.1.4.6), abind(v.1.4-5), compiler(v.4.3.1), bit64(v.4.0.5), withr(v.2.5.0), backports(v.1.4.1), inline(v.0.3.19), shinystan(v.2.6.0), pkgbuild(v.1.4.2), highr(v.0.10), gtools(v.3.9.4), tools(v.4.3.1), httpuv(v.1.6.11), threejs(v.0.3.3), glue(v.1.6.2), callr(v.3.7.3), nlme(v.3.1-162), promises(v.1.2.0.1), checkmate(v.2.2.0), reshape2(v.1.4.4), generics(v.0.1.3), gtable(v.0.3.3), tzdb(v.0.4.0), hms(v.1.1.3), tidygraph(v.1.2.3), xml2(v.1.3.5), utf8(v.1.2.3), pillar(v.1.9.0), ggdist(v.3.3.0), markdown(v.1.7), vroom(v.1.6.3), posterior(v.1.4.1), later(v.1.3.1), lattice(v.0.21-8), bit(v.4.0.5), tidyselect(v.1.2.0), miniUI(v.0.1.1.1), knitr(v.1.43), arrayhelpers(v.1.1-0), gridExtra(v.2.3), V8(v.4.3.3), svglite(v.2.1.1), stats4(v.4.3.1), xfun(v.0.39), bridgesampling(v.1.1-2), matrixStats(v.1.0.0), rstan(v.2.26.22), stringi(v.1.7.12), yaml(v.2.3.7), evaluate(v.0.21), codetools(v.0.2-19), cli(v.3.6.1), RcppParallel(v.5.1.7), shinythemes(v.1.2.0), xtable(v.1.8-4), systemfonts(v.1.0.4), munsell(v.0.5.0), processx(v.3.8.2), jquerylib(v.0.1.4), coda(v.0.19-4), svUnit(v.1.0.6), parallel(v.4.3.1), rstantools(v.2.3.1.1), ellipsis(v.0.3.2), prettyunits(v.1.1.1), dygraphs(v.1.1.1.6), bayesplot(v.1.10.0), Brobdingnag(v.1.2-9), viridisLite(v.0.4.2), mvtnorm(v.1.2-2), scales(v.1.2.1), xts(v.0.13.1), crayon(v.1.5.2), rlang(v.1.1.1), rvest(v.1.0.3) and shinyjs(v.2.1.0)